Stan code for the analysis of my Social and Cultural Dynamics exam
add_predictive_posterior <- function(dataset, model) {
post <- extract.samples(model)
dataset %>%
mutate(theta_ind = purrr::pmap(list(n_ind, group, intensity, vowel), function(n,gro,int,vow) {
iA = post$ia + post$ia_vowel[,vow] + post$ia_group[,gro]
return(inv_logit(iA))
}), theta_gro = purrr::pmap(list(n_gro, group, intensity, vowel), function(n,gro,int,vow) {
gA = post$ga + post$ga_vowel[,vow] + post$ga_group[,gro]
return(inv_logit(gA))
}))
}
my_waic <- function(pred_post) {
pred_post %>%
mutate(
lik_gro = purrr::pmap(list(k_gro, n_gro, theta_gro), function(k,n,theta){
dbinom(k,n,theta)}),
lik_ind = purrr::pmap(list(k_ind, n_ind, theta_ind), function(k,n,theta){
dbinom(k,n,theta)}),
lppd = purrr::map2_dbl(lik_gro, lik_ind, function(lik1, lik2) {
log(mean(c(lik1, lik2)))}),
pwaic = purrr::map2_dbl(lik_gro, lik_ind, ~var2(c(.x, .y)))) %>%
summarise(
WAIC = -2 * (sum(lppd) - sum(pwaic)),
WAIC_SE = sqrt(n() * 2 * var2(-2 * (lppd - pwaic))),
pWAIC = sum(pwaic))
}
dataset = data_frame(
n_gro = full.data$n_gro,
n_ind = full.data$n_ind,
k_gro = full.data$k_gro,
k_ind = full.data$k_ind,
group = full.data$group,
intensity = full.data$intensity,
vowel = full.data$vowel)
do_waic <- function(model, data = dataset) {
add_predictive_posterior(data, model) %>%
my_waic()
}
load("models/intercept.model")
load("models/conf_align_local.model")
load("models/indiscriminate_align_local.model")
load("models/conf_align_global.model")
load("models/cosine.model")
load("models/self_l.model")
load("models/align_l.model")
load("models/synergy_l.model")
load("models/self_entr.model")
load("models/align_entr.model")
load("models/synergy_entr.model")
load("models/self_l_and_local_confidence.model")
load("models/self_l_and_cosine.model")
load("models/self_l_cosine_local_confidence.model")
load("models/cosine_local_confidence.model")
load("models/self_l_x_local_confidence.model")
load("models/self_l_x_cosine.model")
load("models/cosine_x_local_confidence.model")
load("models/self_l_x_cosine_conf_local_full.model")
load("models/self_l_x_cosine_conf_local.model")
load("models/cosine_x_local_confidence_2.model")
load("models/self_l_x_cosine2.model")
load("models/self_l_x_local_confidence2.model")
load("models/self_l_x_local_confidence_only.model")
waic_list <- tribble(
~model, ~waic,
"intercept", do_waic(intercept),
"local confidence alignment", do_waic(conf_align_local),
"local indiscriminate alignment", do_waic(indiscriminate_align_local),
"global confidence alignment", do_waic(conf_align_global),
"word set cosine similarity", do_waic(cosine),
"self L", do_waic(self_l),
"alignment L", do_waic(align_l),
"synergy L", do_waic(synergy_l),
"self entropy", do_waic(self_entr),
"alignment entropy", do_waic(align_entr),
"synergy entropy", do_waic(synergy_entr),
"self L + local confidence", do_waic(self_l_and_local_confidence),
"self L + cosine", do_waic(self_l_and_cosine),
"self L + cosine + local conf", do_waic(self_l_cosine_local_confidence),
"cosine + local confidence", do_waic(cosine_local_confidence),
"self L * conf + cosine", do_waic(self_l_x_local_confidence),
"self L * cosine + local conf", do_waic(self_l_x_cosine),
"cosine * local confidence + self L", do_waic(cosine_x_local_confidence),
"self L * local + self L * cosine", do_waic(self_l_x_local_confidence2),
"self L * cosine + local * cosine", do_waic(self_l_x_cosine2),
"self L * local + cosine * local", do_waic(cosine_x_local_confidence_2),
"self L * cosine * local 2way only",do_waic(self_l_x_cosine_conf_local),
"self L * cosine * local", do_waic(self_l_x_cosine_conf_local_full),
"self L * local", do_waic(self_l_x_local_confidence_only)
) %>% unnest() %>%
arrange(WAIC) %>%
mutate(dWAIC = WAIC - WAIC[1],
weight = exp(-0.5 * dWAIC),
weight = weight / sum(weight)
) %>%
# round the table to 3 digits
mutate_at(vars(-model), function(co) {round(co, 3)})
waic_list
## # A tibble: 24 × 6
## model WAIC WAIC_SE pWAIC dWAIC weight
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 self L * local + self L * cosine 1484.094 112.356 4.906 0.000 0.394
## 2 self L * conf + cosine 1485.220 112.065 4.859 1.126 0.224
## 3 cosine * local confidence + self L 1486.051 111.968 4.842 1.957 0.148
## 4 self L * cosine * local 1486.475 112.913 4.991 2.381 0.120
## 5 self L * cosine + local * cosine 1487.380 112.981 4.890 3.286 0.076
## 6 self L * cosine + local conf 1490.177 112.812 4.847 6.084 0.019
## 7 self L * local + cosine * local 1490.733 112.996 4.894 6.639 0.014
## 8 self L + cosine + local conf 1494.992 113.165 4.766 10.898 0.002
## 9 self L * cosine * local 2way only 1495.961 114.305 4.937 11.867 0.001
## 10 self L * local 1496.377 113.061 4.792 12.283 0.001
## # ... with 14 more rows
# save("waic_list", file="models/waic.Rdata")
pred.post.data <- extract.samples(self_l_x_cosine_conf_local_full) %>%
as.data.frame() %>%
as_data_frame() %>%
select(col_a : col_b7)
pred.post.data %>%
gather() %>%
# group_by()
ggplot(aes(value, fill=key)) +
geom_density(alpha=.4, colour=NA) +
facet_wrap(~ key, scales="free")
pairs(pred.post.data)
# confidence alignment back-scaling
# m = mean(text$local_alignment_confidence)
# s = sd(text$local_alignment_confidence)
#
# data_frame(local_confidence = seq_range(full.data$local_confidence, 10)) %>%
# mutate(samples = rep(list(pred.post.data), n()),
# pred.post = purrr::map2(local_confidence, samples, ~.y$col_a + .y$col_b_confidence * .x),
# local_confidence = (local_confidence * s) + m) %>%
# unnest(pred.post) %>%
# ggplot(aes(local_confidence, pred.post)) +
# geom_line(stat="summary", fun.y=mean) +
# geom_ribbon(stat="summary", fun.data=mean_hpdi, alpha=.3) +
# geom_point(aes(local_confidence, y = .97), data = data.frame(local_confidence = (full.data$local_confidence * s) + m),
# shape="|", size=3) +
# theme_tufte() +
# theme(axis.ticks = element_blank()) +
# labs(x = "Local confidence alignment: Proportion of confidence expressions aligned with previous sentence",
# y = "Predicted collective benefit, mean +- 89% HPDI",
# title = "Counterfactual predictive posterior for explicit metacognition (confidence alignment)",
# subtitle = "Actual dyads' confidence alignment as x-axis") +
# scale_y_continuous(breaks = seq(.96,1.1, .02))
# ggsave("plots/pred_post_local_confidence.png")
posterior = extract.samples(self_l_x_cosine_conf_local_full)
pred <- expand.grid(local_alignment_confidence = seq(-2,2,length.out = 10),
self.l = seq(-2,2, 2),
cosine_word_set = seq(-2,2,2)) %>%
mutate(collective_benefit = purrr::pmap(list(local_alignment_confidence, self.l, cosine_word_set), function(loc, sel, cosi) {
posterior$col_a +
posterior$col_b1 * sel +
posterior$col_b2 * loc +
posterior$col_b3 * cosi +
posterior$col_b4 * loc * sel +
posterior$col_b5 * cosi * sel +
posterior$col_b6 * cosi * loc +
posterior$col_b7 * cosi * loc * sel
})) %>%
plyr::adply(1, function(dat) {
# print(dat)
mean_hpdi(dat$collective_benefit[[1]])
}) %>%
select(-collective_benefit)
pred %>%
ggplot(aes(local_alignment_confidence, ymin=ymin, y=y, ymax=ymax)) +
geom_line() +
geom_ribbon(alpha=.3) +
theme_tufte() +
theme(axis.ticks = element_blank()) +
labs(x = "Local confidence alignment: Proportion of confidence expressions aligned with previous sentence (z-scale)",
y = "Predicted collective benefit, mean +- 89% HPDI",
title = "Counterfactual predictive posterior for the 3 way interaction model",
subtitle = "Facets are z-scaled self L (chat predicability) and cos distance between word sets of interaction partners") +
# scale_y_continuous(limits = c(0,1)) +
coord_cartesian(ylim = c(0,6)) +
facet_grid(cosine_word_set ~ self.l, labeller=label_both)
# ggsave("plots/3way_interaction.png")
posterior = extract.samples(self_l_x_local_confidence2)
pred <- expand.grid(local_alignment_confidence = seq(-2,2,length.out = 10),
self.l = seq(-2,2, 2),
cosine_word_set = seq(-2,2,2)) %>%
mutate(collective_benefit = purrr::pmap(list(local_alignment_confidence, self.l, cosine_word_set), function(loc, sel, cosi) {
posterior$col_a +
posterior$col_b1 * sel +
posterior$col_b2 * loc +
posterior$col_b3 * cosi +
posterior$col_b4 * loc * sel +
posterior$col_b5 * cosi * sel;
})) %>%
plyr::adply(1, function(dat) {
# print(dat)
mean_hpdi(dat$collective_benefit[[1]])
}) %>%
select(-collective_benefit)
pred %>%
ggplot(aes(local_alignment_confidence, ymin=ymin, y=y, ymax=ymax)) +
geom_line() +
geom_ribbon(alpha=.3) +
theme_tufte() +
theme(axis.ticks = element_blank()) +
labs(x = "Local confidence alignment: Proportion of confidence expressions aligned with previous sentence (z-scale)",
y = "Predicted collective benefit, mean +- 89% HPDI",
title = "Counterfactual predictive posterior for: self L * local + self L * cosine",
subtitle = "Facets are z-scaled self L (chat predicability) and cos distance between word sets of interaction partners") +
# scale_y_continuous(limits = c(0,1)) +
coord_cartesian(ylim = c(0,3)) +
facet_grid(self.l ~ cosine_word_set, labeller=label_both)
posterior = extract.samples(self_l_x_local_confidence)
pred <- expand.grid(
local_alignment_confidence = seq(-2,2,length.out = 10),
self.l = seq(-2,2, 2),
# cosine_word_set = seq(-2,2,length.out=10)
# self.l = 0,
# local_alignment_confidence = 0,
cosine_word_set = 0
) %>%
mutate(collective_benefit = purrr::pmap(list(local_alignment_confidence, self.l, cosine_word_set), function(loc, sel, cosi) {
posterior$col_a +
posterior$col_b1 * sel +
posterior$col_b2 * loc +
posterior$col_b3 * cosi +
posterior$col_b4 * loc * sel;
})) %>%
plyr::adply(1, function(dat) {
# print(dat)
mean_hpdi(dat$collective_benefit[[1]])
}) %>%
select(-collective_benefit)
pred %>%
ggplot(aes(local_alignment_confidence, ymin=ymin, y=y, ymax=ymax)) +
geom_line() +
geom_ribbon(alpha=.3) +
theme_tufte() +
theme(axis.ticks = element_blank()) +
labs(x = "Local confidence alignment: Proportion of confidence expressions aligned with previous sentence (z-scale)",
y = "Predicted collective benefit, mean +- 89% HPDI",
title = "Counterfactual predictive posterior for: self L * local + cosine",
subtitle = "Facets are z-scaled self L (chat predicability) and cos distance between word sets of interaction partners") +
# scale_y_continuous(limits = c(0,1)) +
coord_cartesian(ylim = c(0.7,1.4)) +
facet_wrap(~ self.l, labeller=label_both)
# gA = ga + ga_vowel[vowel] + ga_group[group];
# gB = gb + gb_vowel[vowel] + gb_group[group];
# gtheta[i] = gA[i] + gB[i] * intensity[i];
# expand.grid(intensity = unique(full.data$intensity),
# vowel = 1:4,
# group = 1:9) %>%
# mutate(Group = pmap(list(intensity, group, vowel), function(i,g,v){
# post$ga + post$ga_vowel[,v] + post$ga_group[,g] +
# (post$gb + post$gb_vowel[,v] + post$gb_group[,g]) * i
# }), Individual = pmap(list(intensity, group, vowel), function(i,g,v){
# post$ia + post$ia_vowel[,v] + post$ia_group[,g] +
# (post$ib + post$ib_vowel[,v] + post$ib_group[,g]) * i
# })) %>%
# left_join(data.frame(group = 1:9, confidence = full.data$local_confidence)) %>%
# gather("who", "theta", c(Group, Individual)) %>%
# mutate(n = 10,
# theta = purrr::map(theta, ~inv_logit(.))) %>%
# mutate(vowel = c("e", "i", "o", "u")[vowel]) %>%
# unnest(theta) %>%
# ggplot(aes(intensity, theta, group=str_c(group, who), colour=confidence)) +
# # geom_ribbon(aes(group=NA), stat="summary", fun.data=mean_hpdi, alpha=.3) +
# geom_line(stat="summary", fun.y=mean, alpha=.7) +
# geom_line(aes(group=NA), stat="summary", fun.y=mean, colour="red") +
# facet_wrap(~ vowel, labeller = label_both) +
# theme_few() +
# labs(x = "Stimulus intensity (dB)", y = "Proportion correctly identified")